Italian preterms 16S and virome analysis

1 Bacterial taxonomy stats

PSB

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 342 taxa and 61 samples ]
## sample_data() Sample Data:       [ 61 samples by 21 sample variables ]
## tax_table()   Taxonomy Table:    [ 342 taxa by 7 taxonomic ranks ]

total_reads <- sum(sample_sums(PSB))

2 Viral taxonomy stats

total_reads <- sum(sample_sums(PSV))

#Realm-level
# no_realm <- subset_taxa(PSV,Realm=="Unclassified")
# 1-sum(sample_sums(no_realm))/total_reads

#Class
# no_class <- subset_taxa(PSV,Class=="Unclassified")
# 1-sum(sample_sums(no_class))/total_reads

#Order
no_order <- subset_taxa(PSV,Order=="Unclassified")
1-sum(sample_sums(no_order))/total_reads

## [1] 0.0689926

#Family
no_family <- subset_taxa(PSV,Family=="Unclassified")
1-sum(sample_sums(no_family))/total_reads

## [1] 0.005573007

#Genus
no_genus <- subset_taxa(PSV,Genus=="Unclassified")
1-sum(sample_sums(no_genus))/total_reads

## [1] 0

#Species
no_species <- subset_taxa(PSV,Species=="Unclassified")
1-sum(sample_sums(no_species))/total_reads

## [1] 0

#Caudoviricetes
caudovirales <- subset_taxa(PSV,Class=="Caudoviricetes")
sum(sample_sums(caudovirales))/total_reads

## [1] 0.994427

caudovirales_uncl <- subset_taxa(caudovirales,Genus=="Unclassified")
sum(sample_sums(caudovirales_uncl))/sum(sample_sums(caudovirales))

## [1] 1

sum(sample_sums(caudovirales_uncl))/total_reads

## [1] 0.994427

#Classified caudovirales
(sum(sample_sums(caudovirales))-sum(sample_sums(caudovirales_uncl)))/total_reads

## [1] 0

3 Barplots

3.1 Family level

3.2 Bacterial genera plus viral host predictions

4 Abundance heatmaps

4.1 Bacteria

4.1.1 Genus

4.1.2 Species

4.2 Virome

4.2.1 Genus

4.2.2 Species

PSVamp <- ampvis2::amp_load(PSV.no.realm)  

## Phyloseq object provided. Ignoring anything provided for the metadata, taxonomy, fasta, or tree arguments, using those from the phyloseq object instead.

## Warning: Could not find a column named OTU/ASV in otutable, using rownames as
## OTU ID's

## Warning: Could not find a column named OTU/ASV in taxonomy, using rownames as
## OTU ID's

PSVamp$tax[PSVamp$tax == "Unclassified"] <- ""
#rownames(PSVamp$tax) <- gsub("OTU", "vOTU", rownames(PSVamp$tax))

##Virome

amp_vir_species <- amp_heatmap(PSVamp,
            group_by = "Twin_pairs_unique",
            facet_by = "Days",
            plot_values = FALSE,
            tax_show = 30,
            tax_aggregate = "Genus" ,
            tax_add = NULL,
            tax_empty = "best",
            normalise = TRUE,
            color_vector = c("white", "red"),
            plot_colorscale = "sqrt",
            plot_legendbreaks = c(1, 10, 25, 50)
) +
  theme_classic() +
  theme(strip.background = element_rect(colour = "black", fill = "white"),
        strip.text.x = element_text(face = "bold"),
        axis.text.x=element_text(angle = -45, hjust = 0,vjust = 0.2),
        axis.text.y=element_text(face="italic")) +
  scale_y_discrete(labels=function(labels) str_replace_all(labels, c("OTU_" = "unlassified_", "o__" = "", "c__" = "")))



amp_vir_species

4.3 Host virome

5 Upset plot

5.1 Bacteria

5.2 Virome

5.3 Virome host

6 Alpha diversity

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Shannon ~ variable, data = rich)
## 
## $variable
##              diff         lwr       upr     p adj
## 15-1  -0.05465400 -0.88409084 0.7747828 0.9980934
## 30-1   0.04368086 -0.79632248 0.8836842 0.9990587
## 90-1   0.55128566 -0.29487163 1.3974430 0.3208382
## 30-15  0.09833487 -0.48059748 0.6772672 0.9694574
## 90-15  0.60593967  0.01811383 1.1937655 0.0409439
## 90-30  0.50760480 -0.09503884 1.1102484 0.1276992

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Shannon ~ variable, data = rich)
## 
## $variable
##              diff        lwr       upr     p adj
## 15-1  -0.28313685 -0.8951040 0.3288303 0.6133665
## 30-1  -0.01146464 -0.6422659 0.6193367 0.9999594
## 90-1  -0.39036878 -1.0594346 0.2786970 0.4177807
## 30-15  0.27167221 -0.3402949 0.8836394 0.6443038
## 90-15 -0.10723193 -0.7585710 0.5441071 0.9719606
## 90-30 -0.37890414 -1.0479699 0.2901617 0.4442411

6.1 Bacterial and viral shannon diversity

7 Beta diversity

7.1 Days

8 Within vs between twin pairs bray curtis dissimilarity

8.1 Individual time points

8.1.1 Day 1

8.1.2 Day 15

8.1.3 Day 30

8.1.4 Day 90

9 Effect size

9.1 Bacteria

## [1] 137

9.1.1 Db anova - constrained by time point

9.2 Virome

## [1] 576

9.2.1 Db anova - constrained for time point

9.3 Host virome

## [1] 305

9.3.1 Db anova - constrained for time point

9.4 Host virome - Genus

## [1] 305

##  [1] "Twin_pairs"                         "Days"                              
##  [3] "Gestational_age"                    "Birth_weight"                      
##  [5] "Small_for_gestational_age"          "Antenatal_steroids"                
##  [7] "Delivery_mode"                      "Intrapartum_antibiotic_prophylaxis"
##  [9] "Late_onset_sepsis"                  "Necrotizing_enterocolitis"         
## [11] "Patent_ductus_arteriosus"           "Bronchopulmonary_dysplasia"        
## [13] "Antibiotics_at_birth"               "Total_days_antibiotics_therapy"

9.4.1 Db anova - constrained for time point

9.5 Combined effect size plots

9.5.1 With viral hosts

9.5.2 Bacteria vs vira only

10 Rabuplot

10.1 Time

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

PSV.host@sam_data$Days <- factor(PSV.host@sam_data$Days, levels = c("90","60","30","15","1"))
PSB@sam_data$Days <- factor(PSB@sam_data$Days, levels = c("90","60","30","15","1"))

rabu.vir <- rabuplot(
  PSV.host,
  violin = TRUE,
  predictor = "Days",
  p_adjust = TRUE,
  N_taxa = 10,
  colors = col_fil[1:4],
  italic_names = FALSE) +
  ggtitle("Viral hosts")

labels <- levels(factor(rabu.vir$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.vir <- rabu.vir + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())

rabu.bac <- rabuplot(
  PSB,
  violin = TRUE,
  predictor = "Days",
  p_adjust = TRUE,
  N_taxa = 10,
  colors = col_fil[1:4],
  italic_names = FALSE) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria")

labels <- levels(factor(rabu.bac$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.bac <- rabu.bac + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())

ggarrange(rabu.bac,rabu.vir, common.legend = TRUE, legend = "bottom")

10.2 Bronchopulmonary_dysplasia

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.bpd <- rabuplot(
  PSV.host,
  violin = TRUE,
  predictor = "Bronchopulmonary_dysplasia",
  legend_title = "Bronchopulmonary dysplasia",
  p_adjust = TRUE,
  #p_stars = TRUE,
  N_taxa = 10,
  colors = col_fil[1:2],
  italic_names = FALSE) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

labels <- levels(factor(rabu.vir.bpd$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.vir.bpd <- rabu.vir.bpd + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())


rabu.bac.bpd <- rabuplot(
  PSB,
  violin = TRUE,
  predictor = "Bronchopulmonary_dysplasia",
  legend_title = "Bronchopulmonary dysplasia",
  p_adjust = TRUE,
  #p_stars = TRUE,
  N_taxa = 10,
  colors = col_fil[1:2],
  italic_names = FALSE) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

labels <- levels(factor(rabu.bac.bpd$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.bac.bpd <- rabu.bac.bpd + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())
  

ggarrange(rabu.bac.bpd,rabu.vir.bpd, common.legend = TRUE, legend = "bottom")

10.3 Intrapartum antibiotics

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.tdab <- rabuplot(
  PSV.host,
  violin = TRUE,
  predictor = "Intrapartum_antibiotic_prophylaxis",
  legend_title = "Intrapartum antibiotic prophylaxis",
  p_adjust = TRUE,
  N_taxa = 10,
  colors = col_fil[1:3],
  italic_names = FALSE) +
  #theme(legend.position = "none") +
  ggplot2::ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

labels <- levels(factor(rabu.vir.tdab$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.vir.tdab <- rabu.vir.tdab + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())

rabu.bac.tdab <- rabuplot(
  PSB,
  violin = TRUE,
  predictor = "Intrapartum_antibiotic_prophylaxis",
  legend_title = "Intrapartum antibiotic prophylaxis",
  p_adjust = TRUE,
  N_taxa = 10,
  colors = col_fil[1:3],
  italic_names = FALSE) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

labels <- levels(factor(rabu.bac.tdab$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.bac.tdab <- rabu.bac.tdab + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())

ggarrange(rabu.bac.tdab,rabu.vir.tdab, common.legend = TRUE, legend = "bottom")

10.4 Total days antibiotics treatment

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.iap <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Total_days_antibiotics_therapy", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.iap <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Total_days_antibiotics_therapy", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.iap,rabu.vir.iap, common.legend = TRUE, legend = "bottom")

10.5 Antibiotics at birth

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.tdab <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Antibiotics_at_birth", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.tdab <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Antibiotics_at_birth", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.tdab,rabu.vir.tdab, common.legend = TRUE, legend = "bottom")

10.6 Small for gestational age

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.sga <- rabuplot(
  PSV.host,
  violin = TRUE,
  predictor = "Small_for_gestational_age",
  legend_title = "Small for gestational age",
  p_adjust = TRUE,
  #p_stars = TRUE,
  N_taxa = 10,
  colors = col_fil[3:4],
  italic_names = FALSE) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

labels <- levels(factor(rabu.vir.sga$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.vir.sga <- rabu.vir.sga + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())

rabu.bac.sga <- rabuplot(
  PSB,
  violin = TRUE,
  predictor = "Small_for_gestational_age",
  legend_title = "Small for gestational age",
  p_adjust = TRUE,
  #p_stars = TRUE,
  N_taxa = 10,
  colors = col_fil[3:4],
  italic_names = FALSE) +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

labels <- levels(factor(rabu.bac.sga$data$variable))
labels[labels != "Other genus"] <- paste0("*", labels[labels != "Other genus"], "*")

rabu.bac.sga <- rabu.bac.sga + scale_x_discrete(labels = labels) + theme(axis.text.y = ggtext::element_markdown())

ggarrange(rabu.bac.sga,rabu.vir.sga, common.legend = TRUE, legend = "bottom")

11 RCCA - trans kingdom interactions

11.1 Conversion and screeplot

#Subset samples present in both
vir_ids <- PSV@sam_data$description %>% str_split_i("_",-1)
PSV@sam_data$id <- vir_ids

PSB.CSS.both <- subset_samples(PSB,description %in% vir_ids)
#PSV.CSS.both <- subset_samples(PSV.CSS,vir_ids %in% PSB.CSS.both@sam_data$description)
PSV.CSS.host.both <- subset_samples(PSV.CSS.host,vir_ids %in% PSB.CSS.both@sam_data$description)

##Transpose and convert to dataframes
botu <- data.frame(PSB.CSS.both@otu_table)
#Select last non-empty taxonomic rank
tax.bac <- as.data.frame(PSB.CSS.both@tax_table) %>% dplyr::mutate_all(~na_if(., '')) %>% apply(.,1, function(x) tail(stats::na.omit(x), 1))
rownames(botu) <- paste(rownames(botu),tax.bac, sep = "-")

votu <- data.frame(PSV.CSS.host.both@otu_table)
#Select last non-empty taxonomic rank
tax.vir <- as.data.frame(PSV.CSS.host.both@tax_table) %>% dplyr::mutate_all(~na_if(., '')) %>% apply(.,1, function(x) tail(stats::na.omit(x), 1))
rownames(votu) <- paste(rownames(votu),tax.vir, sep = "-")


  ###########Convert to matrix and transpose - Maybe skip?
t(as.matrix.data.frame(botu))->tmbacteria
  #Save OTU Taxonomy from header
tmcolnames <- colnames(tmbacteria)
  #Remove first line
rownames(tmbacteria) <- NULL
  #Apply as.numeric to matrix(will remove header as it is in characters)
nbacteria <- mapply(tmbacteria, FUN=as.numeric)
  #Convert back into matrix and add colnames
mbacteria <- matrix(nbacteria, ncol=NCOL(tmbacteria), nrow=NROW(tmbacteria))
colnames(mbacteria) <- tmcolnames

  
    ##Repeat for viral OTU table
t(as.matrix.data.frame(votu))->tmvirome
tmcolnames <- colnames(tmvirome)
rownames(tmvirome) <- NULL
nvirome <- mapply(tmvirome, FUN=as.numeric)
mvirome <- matrix(nvirome, ncol=NCOL(tmvirome), nrow=NROW(tmvirome))
colnames(mvirome) <- tmcolnames
    ##Alternatively use: data.matrix(frame, rownames.force = NA)

    ##Converting and removing near-zero variance values
bac <- mbacteria[,-c(1)]
nzv <- nearZeroVar(bac)
bac <- bac[, -nzv$Position]
           
           
vir <- mvirome[,-c(1)]
#nzv <- nearZeroVar(vir) # If any
#vir <- vir[, -nzv$Position]


  #PCA
library(vegan)
mod_bac <- rda(bac, scale = TRUE)
plot(mod_bac, scaling = 3)

screeplot(mod_bac)

mod_vir <- rda(vir, scale = TRUE)
plot(mod_vir, scaling = 3)

screeplot(mod_vir)

scaling = "shrinkage"
  #Apply scaling 
result_1comp <- rcc(bac, vir, method = scaling, ncomp = 1)
result_2comp <- rcc(bac, vir, method = scaling, ncomp = 2)
result_3comp <- rcc(bac, vir, method = scaling, ncomp = 3)
result_4comp <- rcc(bac, vir, method = scaling, ncomp = 4)

ggarrange(plot(result_1comp, type = "barplot"),
          plot(result_2comp, type = "barplot"),
          plot(result_3comp, type = "barplot"),
          plot(result_4comp, type = "barplot")
          ,nrow = 2)

## $`1`

## 
## $`2`

## 
## attr(,"class")
## [1] "list"      "ggarrange"

# produce a heat map of the cross correlation matrix
imgCor(bac, vir, sideColors = c("purple", "green")) 

11.2 Model tuning

# set grid search values for each regularisation parameter
grid1 <- seq(0.001, 0.2, length = 20) 
grid2 <- seq(0.001, 0.2, length = 20)

# optimise the regularisation parameter values
sl(cv.tune.rcc, tune.rcc(bac, vir, grid1 = grid1, grid2 = grid2, 
                                   validation = "loo"),overwrite = F)

# examine the results of CV tuning
plot(cv.tune.rcc)

cv.tune.rcc

## 
## Call:
##  tune.rcc(X = bac, Y = vir, grid1 = grid1, grid2 = grid2, validation = "loo") 
## 
##   lambda1 =  0.2 , see object$opt.lambda1
##   lambda2 =  0.001 ,  see object$opt.lambda2
##  CV-score =  0.2734331 , see object$opt.score

#manually set lambda scores
#opt.l1 <- 1.7
#opt.l2 <- 0.003
opt.l1 <- cv.tune.rcc$opt.lambda1
opt.l2 <- cv.tune.rcc$opt.lambda2

11.3 Model and component selection

shrink.rcc <- rcc(bac, vir, method = "shrinkage")

CV.rcc <- rcc(bac, vir, method = "ridge", 
                         lambda1 = opt.l1,
                         lambda2 = opt.l2) 

# barplot of cross validation method rCCA canonical correlations
plot(CV.rcc, type = "barplot", main = "Cross Validation") 

# barplot of shrinkage method rCCA canonical correlations
plot(shrink.rcc, type = "barplot", main = "Shrinkage")

11.4 Run model

#Perfrom optimised CV rCCA

result.CV <- rcc(bac, vir, method = "ridge",
                         lambda1 = opt.l1,
                         lambda2 = opt.l2,
              ncomp=3)

result.raw <- rcc(bac, vir, method = "shrinkage",ncomp=3)

#Remove OTU IDs to only show taxonomy
result <- result.raw

result$names$colnames$X <- gsub(".*\\-", "", result$names$colnames$X) %>% make.unique(sep="_")
result$names$colnames$Y <- gsub(".*\\-", "", result$names$colnames$Y) %>% make.unique(sep="_")


# colnames(result$X) <- gsub(".*\\-", "", colnames(result$X))
# colnames(result$Y) <- gsub(".*\\-", "", colnames(result$Y))

# rownames(result$loadings$X) <- gsub(".*\\-", "", rownames(result$loadings$X))
# rownames(result$loadings$Y) <- gsub(".*\\-", "", rownames(result$loadings$Y))

11.5 Visualization

11.5.1 Heatmaps

#Defne pallette to only add colour to only high
##Generate new pallette
palette <- hcl.colors(200, palette = "Zissou 1", alpha = NULL, rev = FALSE, fixup = TRUE) #Seclect palette https://developer.r-project.org/Blog/public/2019/04/01/hcl-based-color-palettes-in-grdevices/
palette[61:140] <- "white" #Select which part of pallette to replace with white

# high4_2c <- cim(result, comp = 1:2, ylab = "Bacteria", xlab = "Virome", margins = c(12,18),cutoff = 0.4, color = palette)
# high45_2c <- cim(result, comp = 1:2, ylab = "Bacteria", xlab = "Virome", margins = c(12,18),cutoff = 0.45, color = palette)
#high5_3c <- cim(result, comp = 1:2, ylab = "Bacteria", xlab = "Virome", margins = c(12,18),cutoff = 0.5, color = palette)

#pdf(filename="high5_3c.pdf", res = 600, width = 4000, height = 4000)
#high5_3c
#dev.off()

######

high6_3c <- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(7,9),cutoff = 0.6, color = palette)

high5_3c <- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(7,9),cutoff = 0.5, color = palette)

high45_3c<- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(15,20),cutoff = 0.45, color = palette)

high4_3c <- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(15,20),cutoff = 0.4, color = palette)

#mat5 <- high5_2c$mat
mat6_3c <- high6_3c$mat
mat45_3c <- high45_3c$mat
mat4_3c <- high4_3c$mat
  
 #Output Matrix
write.table(mat45_3c, "Matrix_high_3-Component_0.45.txt", sep="\t")

#Make linear plots of all bac and vir 
bacmat6_3c <- bac[,rownames(mat6_3c)]
virmat6_3c <- vir[,colnames(mat6_3c)]
bacvirmat6_3c <- cbind(bacmat6_3c,virmat6_3c)
bac45_3c <- bac[,rownames(mat45_3c)]
vir45_3c <- vir[,colnames(mat45_3c)]
bacvirmat45_3c <- cbind(bac45_3c,vir45_3c)

#ggpairs
library(GGally)
#ggpairs(bacvir)
#ggpairs(bacvir45,columns=3:13)
#ggpairs(as.data.frame(bacvir45))
sl(pairs45,ggpairs(as.data.frame(bacvirmat45_3c)))
sl(pairs6,ggpairs(as.data.frame(bacvirmat6_3c)))

# pairs45 +
#   theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) 

#pairs6 +
#   theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) 

11.5.2 Network plots

11.5.2.1 Shrinkage

# Reformat for gephi

## Save blocks as polygon to determine shape in network

result$names$Polygon <- result$names$data
result$names$Polygon <- replace(result$names$data, c("bac","vir"),c(3,4)) %>% as.integer()

## Remove ubique markers to allow shared coloring in Gephi
result$names_simple$colnames$X <- gsub("*.\\_", "", result$names$colnames$X)
result$names_simple$colnames$Y <- gsub("*.\\_", "", result$names$colnames$Y)

#Show network

network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cex.node.name = 14,
        cutoff = 0.3,
        color.node = col_fil[1:2]) 

network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.4,
        color.node = col_fil[1:2]) 

network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.5,
        color.node = col_fil[1:2]) 

#Export network result
network.res <- network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.30,
        color.node = col_fil[1:2])

library(igraph)

#network.res 

write.graph(network.res$gR, file = "network.shrinkage0.3.gml", format = "gml")

11.5.2.2 CV

result.CV$names_simple$colnames$X <- gsub("*.\\_", "", result.CV$names$colnames$X)
result.CV$names_simple$colnames$Y <- gsub("*.\\_", "", result.CV$names$colnames$Y)

#Show network
network(result.CV, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.3,
        color.node = col_fil[1:2]) 

network(result.CV, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.35,
        color.node = col_fil[1:2]) 

# network(result.CV, comp = 1:3, interactive = FALSE,
#         lwd.edge = 2,
#         cutoff = 0.5,
#         color.node = col_fil[1:2]) 


library(igraph)
network.res <- network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.30,
        color.node = col_fil[1:2])

write.graph(network.res$gR, file = "network.CV0.3.gml", format = "gml")

12 Figures - arrange and save

12.1 Figure 1

12.2 Figure S1

12.3 Figure 2

12.4 Figure S2

12.5 Figure S3

12.6 Figure 3

12.7 Figure S4

13 Session info

sessionInfo()

## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Danish_Denmark.utf8  LC_CTYPE=Danish_Denmark.utf8   
## [3] LC_MONETARY=Danish_Denmark.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=Danish_Denmark.utf8    
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] igraph_1.6.0                GGally_2.2.0               
##  [3] rabuplot_0.0.1.09           UpSetR_1.4.0               
##  [5] VennDiagram_1.7.3           futile.logger_1.4.3        
##  [7] reshape2_1.4.4              metagenomeSeq_1.38.0       
##  [9] RColorBrewer_1.1-3          glmnet_4.1-8               
## [11] Matrix_1.6-5                limma_3.52.4               
## [13] readxl_1.4.3                phyloseq_1.40.0            
## [15] mixOmics_6.20.0             MASS_7.3-60.0.1            
## [17] usdm_2.1-7                  terra_1.7-65               
## [19] directlabels_2024.1.21      vegan_2.6-4                
## [21] lattice_0.22-5              permute_0.9-7              
## [23] viridis_0.6.5               viridisLite_0.4.2          
## [25] DESeq2_1.36.0               SummarizedExperiment_1.26.1
## [27] Biobase_2.56.0              MatrixGenerics_1.8.1       
## [29] matrixStats_1.2.0           GenomicRanges_1.48.0       
## [31] GenomeInfoDb_1.32.4         IRanges_2.30.1             
## [33] S4Vectors_0.34.0            BiocGenerics_0.42.0        
## [35] data.table_1.14.10          RVAideMemoire_0.9-83-7     
## [37] cowplot_1.1.3               plyr_1.8.9                 
## [39] ampvis2_2.8.6               rstatix_0.7.2              
## [41] ggsci_3.0.0                 ggpubr_0.6.0               
## [43] lubridate_1.9.3             forcats_1.0.0              
## [45] stringr_1.5.1               dplyr_1.1.4                
## [47] purrr_1.0.2                 readr_2.1.5                
## [49] tidyr_1.3.1                 tibble_3.2.1               
## [51] ggplot2_3.5.1               tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3          ragg_1.2.7              bit64_4.0.5            
##   [4] knitr_1.45              multcomp_1.4-25         DelayedArray_0.22.0    
##   [7] KEGGREST_1.36.3         MicrobiotaProcess_1.8.2 RCurl_1.98-1.14        
##  [10] generics_0.1.3          lambda.r_1.2.4          TH.data_1.1-2          
##  [13] usethis_2.2.2           RSQLite_2.3.5           commonmark_1.9.1       
##  [16] bit_4.0.5               tzdb_0.4.0              ggtreeExtra_1.6.1      
##  [19] xml2_1.3.6              httpuv_1.6.14           xfun_0.41              
##  [22] hms_1.1.3               jquerylib_0.1.4         evaluate_0.23          
##  [25] promises_1.2.1          DEoptimR_1.1-3          fansi_1.0.6            
##  [28] caTools_1.18.2          DBI_1.2.1               geneplotter_1.74.0     
##  [31] htmlwidgets_1.6.4       rARPACK_0.11-0          ellipsis_0.3.2         
##  [34] RSpectra_0.16-1         ggnewscale_0.4.9        backports_1.4.1        
##  [37] markdown_1.12           annotate_1.74.0         libcoin_1.0-10         
##  [40] vctrs_0.6.5             remotes_2.4.2.1         here_1.0.1             
##  [43] abind_1.4-5             httr2_1.0.0             cachem_1.0.8           
##  [46] withr_3.0.0             robustbase_0.99-2       treeio_1.20.2          
##  [49] getopt_1.20.4           cluster_2.1.6           ape_5.7-1              
##  [52] lazyeval_0.2.2          crayon_1.5.2            ellipse_0.5.0          
##  [55] genefilter_1.78.0       pkgconfig_2.0.3         labeling_0.4.3         
##  [58] nlme_3.1-164            pkgload_1.3.4           gitcreds_0.1.2         
##  [61] devtools_2.4.5          rlang_1.1.3             lifecycle_1.0.4        
##  [64] miniUI_0.1.1.1          sandwich_3.1-0          cellranger_1.1.0       
##  [67] rprojroot_2.0.4         aplot_0.2.2             raster_3.6-26          
##  [70] lpsymphony_1.24.0       carData_3.0-5           boot_1.3-28.1          
##  [73] Rhdf5lib_1.18.2         zoo_1.8-12              png_0.1-8              
##  [76] bitops_1.0-7            KernSmooth_2.23-22      rhdf5filters_1.8.0     
##  [79] Biostrings_2.64.1       blob_1.2.4              shape_1.4.6            
##  [82] coin_1.4-3              gridGraphics_0.5-1      ggsignif_0.6.4         
##  [85] scales_1.3.0            memoise_2.0.1           magrittr_2.0.3         
##  [88] gplots_3.1.3            zlibbioc_1.42.0         compiler_4.2.1         
##  [91] lme4_1.1-35.3           ggstar_1.0.4            cli_3.6.2              
##  [94] ade4_1.7-22             XVector_0.36.0          urlchecker_1.0.1       
##  [97] lmerTest_3.1-3          patchwork_1.2.0         formatR_1.14           
## [100] mgcv_1.9-1              tidyselect_1.2.1        stringi_1.8.4          
## [103] textshaping_0.3.7       highr_0.10              yaml_2.3.8             
## [106] locfit_1.5-9.8          ggrepel_0.9.5           sass_0.4.8             
## [109] ggstats_0.5.1           tools_4.2.1             timechange_0.3.0       
## [112] parallel_4.2.1          rstudioapi_0.15.0       foreach_1.5.2          
## [115] optparse_1.7.4          gridExtra_2.3           farver_2.1.2           
## [118] digest_0.6.35           shiny_1.8.0             ggtext_0.1.2           
## [121] quadprog_1.5-8          Rcpp_1.0.12             gridtext_0.1.5         
## [124] car_3.1-2               broom_1.0.5             later_1.3.2            
## [127] httr_1.4.7              AnnotationDbi_1.58.0    Wrench_1.14.0          
## [130] colorspace_2.1-0        XML_3.99-0.16.1         fs_1.6.4               
## [133] splines_4.2.1           yulab.utils_0.1.4       tidytree_0.4.6         
## [136] gh_1.4.0                sp_2.1-3                multtest_2.52.0        
## [139] ggplotify_0.1.2         systemfonts_1.0.5       plotly_4.10.4          
## [142] sessioninfo_1.2.2       xtable_1.8-4            nloptr_2.0.3           
## [145] futile.options_1.0.1    jsonlite_1.8.8          ggtree_3.4.4           
## [148] corpcor_1.6.10          modeltools_0.2-23       ggfun_0.1.4            
## [151] R6_2.5.1                profvis_0.3.8           metagMisc_0.0.4        
## [154] pillar_1.9.0            htmltools_0.5.7         mime_0.12              
## [157] minqa_1.2.7             Maaslin2_1.10.0         glue_1.7.0             
## [160] fastmap_1.1.1           BiocParallel_1.30.4     codetools_0.2-19       
## [163] pcaPP_2.0-4             pkgbuild_1.4.4          mvtnorm_1.2-4          
## [166] utf8_1.2.4              bslib_0.6.1             numDeriv_2016.8-1.1    
## [169] curl_5.2.0              gtools_3.9.5            biglm_0.9-2.1          
## [172] survival_3.5-7          rmarkdown_2.25          biomformat_1.24.0      
## [175] munsell_0.5.1           rhdf5_2.40.0            GenomeInfoDbData_1.2.8 
## [178] iterators_1.0.14        gtable_0.3.5